Free energy functionals for efficient phase field crystal modeling of structural phase transformations 
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The phase field crystal (PFC) method has emerged as a promising technique for modeling materials with 
atomistic resolution on mesoscopic time scales. The approach is numerically much more efficient than classical 
density functional theory (CDFT), but its single mode free energy functional only leads to lattices with triangular 
(2D) or BCC (3D) symmetries. By returning to a closer approximation of the CDFT free energy functional, we 
develop a systematic construction of two-particle direct correlation functions that allow the study of a broad 
class of crystalline structures. This construction examines planar spacings, lattice symmetries, planar atomic 
densities and the atomic vibrational amplitude in the unit cell of the lattice and also provides control parameters 
for temperature and anisotropic surface energies. The power of this new approach is demonstrated by two 
examples of structural phase transformations. 
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Solid-state transformations form the basis of many prob- 
lems in materials science and physical metallurgy Jl]]. They 
involve complex structural changes between parent and 
daughter phases and couple atomic-scale elastic and plas- 
tic effects with diffusional processes. These phenomena are 
presently impossible to compute at experimentally relevant 
time scales using molecular dynamics simulations. On the 
other hand, meso-scale continuum models wash out most of 
the relevant atomic scale physics that leads to elasticity, plas- 
ticity, defect interactions and grain boundary nucleation and 
migration. Phase field studies of precipitate and ledge growth 
12j-|4|] must thus re-introduce these effects phenomenologi- 
cally. There is presently no continuum method to efficiently 
simulate solid-state transformations on diffusional time scales 
that self-consistently computes elastic and plastic effects at 
the atomic scale. 

Classical density function theory (CDFT) provides a for- 
malism that can accurately describe the emergence of crys- 
talline order from a liquid or solid phase through a coarse- 
grained density field |5|]. Unfortunately, this approach re- 
quires very high spatial resolution and is highly inefficient 
for dynamical calculations J(J. A recent model, coined the 
phase field crystal (PFC) model, has been gaining widespread 
recognition as a hybrid method between CDFT and traditional 
phase field methods. PFC models capture most of the essential 
physics of CDFT without having to resolve atomically sharp 
atomic density peaks J7-1J.]. In spite of their successes, how- 
ever, only PFC free energy functionals minimized by trian- 
gular (2D) or BCC (3D) lattices have been seriously studied. 
While there have been attempts to extend the PFC model to 
describe other crystal symmetries such as square lll2[ll3ll and 
FCC SEE] lattices, these have been somewhat ad-hod and 
not self-consistently connected to material properties. PFC 
modeling is presently lacking a generalized free energy for- 
mulation that allows the study of important phase transforma- 
tions between different crystalline states. 

This Letter proposes a new CDFT/PFC model based on 



two-point direct correlation function that is systematically 
constructed to be minimized by arbitrary crystalline states. 
The utility and efficiency of our density functional model is 
demonstrated by dynamically simulating the growth of solid 
phases into a liquid and the nucleation of precipitate phases 
within a parent phase of different crystallographic symmetry. 
We begin with the CDFT free energy functional, 



AF[n(r)} = AF ld [n(r)} + AF ex [n(f)} 



(1) 



where AF is the free energy difference with respect to a ref- 
erence state at density po 15|], and has contributions from a 
non-interacting term AFid and an excess energy AF ex . The 
latter is responsible for the formation of structured phases. 

Eq. ([U is expanded around the reference density using the 
dimensionless number density n(f) = p(r)/p — 1, where 
p(f) is the local density. The non-interacting contribution is 
approximated by expanding the Helmholtz free energy of an 
ideal gas, 



AF. ld = Po k B T I dr(l + n(f))ln(l + n(r)) - n{r) 
n(r) 2 n(f) 3 n(r) 4 " 1 
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The excess energy is expanded as in CDFT |3l and earlier PFC 
models ^[lfj up to second order correlations, 



AF, 



p k B T 



drn{r) / df (C 2 (|r - r'\)n(?))(?) 



where C2(\r — r'\) is the two particle direct correlation func- 
tion. The dimensionless number density field n(r) is evolved 
in time using conserved dissipative dynamics, 



dn(r) 
dt 



SAF 

8n(r) 
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where M is a kinetic mobility parameter and rj is conserved 
random noise. For increased numerical efficiency, Eq. (O is 
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FIG. 1. (a) Direct pair correlation function for a square lattice at 
two temperatures (a — 0.82, solid and a = 0.5, dotted, see text). 
The two peaks in the correlation function each represent one of two 
planes in the unit cell with interplanar spacings Ai and A2 (see in- 
set), (b) Correlation function for a simple cubic (SC) lattice at two 
different temperatures. In 3D, a cubic lattice is represented by three 
planes in the unit cell. The simple cubic cell has two planar spacings 
in common with a BCC lattice and requires the third peak to stabilize 
the structure. 



integrated in reciprocal space using a semi-implicit technique 

flu! ■ 

Periodic structures emerge in the phase field crystal model 
through the input of a direct correlation function into the ex- 
cess free energy. Here we shall construct correlation functions 
to create free energy minima for arbitrary crystalline lattices 
by expanding the free energy functional about a solid state 
reference density 111 711 . The original PFC model uses a single 
mode correlation function; our correlation functions include 
multiple modes, each of which corresponds to the spacing be- 
tween atomic planes within a particular unit cell. Figure Q] 
shows two examples a the direct pair correlation function in 
reciprocal space that generate stable square (a) and simple cu- 
bic (SC) (b) lattices when the corresponding free energy is 
minimized. The functions Cz(k) are constructed in recipro- 
cal space by combining multiple peaks whose position, am- 
plitudes and width are determined as follows: 

The number of peaks in the correlation function and their 
reciprocal space positions are determined by the unit cell of 
the desired lattice. Every plane within the unit cell produces 
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FIG. 2. Effect of peak width on in the correlation functions on planar 
interface width W%. Left: Three values of on and their resulting 
interface widths. Top Right: Peak shapes corresponding the interface 
widths to the left. Bottom Right: Dependence of simulated interface 
width on the fourier space correlation peak width (see text). 



a family of peaks derived from the inter-planar spacings. For 
a 2D square lattice, for instance, the unit cell contains two 
families of planes, {10} and {11}, with spacings Ai and A2 as 
shown in the inset of Figure QJ a). For a perfect crystal, there 
are an infinite number of peaks located at integer multiples 
of the wavevector fc, = 2ir/\i, where i denotes the plane. 
For square and cubic lattices, it is sufficient to keep only the 
lowest frequency mode for each family of these peaks in the 
correlation function for each plane in the unit cell. 

Temperature dependence enters the correlation function via 
modulation of the peak heights through a Debye-Waller factor, 
denoted Ej,. As thermal motion causes atoms to vibrate about 
their lattice positions, the reciprocal space lattice amplitudes 
fall off by e~ CT where a is related to the vibrational am- 
plitude of the atoms and thus temperature. The peak height 
is additionally influenced by the (dimensionless) atomic den- 
sity pi within a plane and the number of planes /?j in each 
family. This effect can be seen when integrating out the an- 
gular dependence in the reciprocal lattice. Here we modify 
the Debye-Waller Factor to include both of these effects as 

Pi = exp [~ 2 .i. ]- For example, in a square lattice the fam- 
ilies of {11} and {10} planes each consist of 4 sets of planes 
(ie. the {11} family contains (11), (11), (11), and (11)) and 
therefore /3n = /3io = 4. The (ll)-plane has an atomic den- 
sity pn = and the (10)-plane has a density of pxo = 1- 
The peaks in the reciprocal space correlation function are 
represented by Gaussians cxp|— k 2a i' \ of finite width ai 
rather than the 5— peaks of a perfect lattice. The parameter 
ai provides a way to account for changes in the free energy 
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FIG. 3. (a) Phase diagram for a square lattice correlation function 
showing coexistence between liquid, solid and triangular phases. 
The inset contains stripes of square and triangular lattices in a liquid 
phase quenched to a temperature a = 0.79. The square phase den- 
sity was initialized to n sq = 0.067 and the triangular phase density 
to ntri = 0.076. The mean density of the system is set to n = 0.07 
and the surface energy parameters are a.\ = y2 and 012 = 1. (b) 
Phase diagram for a FCC correlation function. The insets show the 
peritectic point and the energy curves for the liquid, FCC and BCC 
states at temperature of a = 0.06. 

due to interfaces, defects and strain. The effect of varying a, t 
on the width of a liquid-solid interface is illustrated in Figure 
12 Increasing the value of ct<j leads to a decrease in the inter- 
face width, which directly affects the surface energy as shown 
previously by Majaniemi and Provatas lfl8ll . In fact, the rela- 
tionship between the interface width and the peak width in the 
correlation function is 1 jai oc W, as illustrated in Figure |2 
This result can be arrived at by the inverse fourier transform of 
C2(k) and agrees well with the traditional phase field models, 
which incorporate surface energy through square gradients in 
the order parameter. Note also that oti incorporates temper- 
ature dependence in the surface energy through the Debye- 
Waller shift of the peak height. 

In summary, each family of planes i in the unit cell con- 
tributes a peak to the direct pair correlation function of the 
form 

a 2 k 2 _(k-k i ) 2 

C 2 {k) l = e-^e 2a 't . (5) 



The correlation functions seen in Figure Q] are the envelope 
of the superposition of all relevant peaks for the crystal struc- 
ture. This construction described above provides a a general 
and robust method by which to generate desired crystal struc- 
tures, which incorporate temperature dependence through the 
Debye- Waller factor. This technique also allows for simple 
control over the surface energy, which can be independently 
tuned to incorporate anisotropy between planar facets. 

By substituting the correlation functiond into the free en- 
ergy Eq. (Q~|), the phase diagram is obtained at different temper- 
atures, using the free energy curves obtained for each phase. 
The energy curves are calculated by an iterative relaxation 
technique for each structure using the minimal kernel for the 
structure of interest. For the liquid state, the energy is calcu- 
lated by imposing a constant density field, and the energy per 
unit volume is calculated by numerically integrating Eq. (Q]i. 
For the solid state, the structures corresponding to the corre- 
lation function are fit using a Gaussian density field. For ex- 
ample, for the correlation function of Fig. 02a), density fields 
for square and triangular lattices with lattice spacings of 1 and 
l/v^3 are constructed and relaxed using Eq. to allow the 
density peaks to obtain an amplitude corresponding to the lo- 
cal minimum energy state. This process is repeated for a series 
of mean densities to produce energy-density curves for each 
phase at a given quench. An example of these curves is shown 
as an inset to Fig.|3b) for an FCC free energy functional. The 
double tangent construction is performed for many quenches 
giving the phase diagrams in Fig. a) for a square correlation 
function and in Fig. 0b) for an FCC correlation function. 

In both 2D and 3D, we find that the phase that possesses 
the symmetry of the underlying correlation function emerges 
at low temperatures and intermediate densities. As the temper- 
ature is increased, the lowest frequency mode becomes dom- 
inant and for high enough densities, the square (FCC) crystal 
transforms into a triangular (BCC) crystal. For lower densi- 
ties, the square (FCC) crystal melts into a homogeneous liq- 
uid. Note that our model also predicts a peritetic point in both 
dimensions, where two solid phases coexist with a liquid. This 
remarkable feature opens a window into the study of structural 
phase transformations and structural phase coexistence. Co- 
existence between FCC crystals with the liquid phase is also 
easily obtained. Similarly, the correlation function shown in 
Fig.[TJb) produces coexistence of a liquid and SC phase. 

Our new model is applied here to two important exam- 
ples of nonequilibrium solid state processese. The first is the 
growth of two structurally different lamellae (stripes) into a 
liquid phase as illustrated in the inset of Fig.[5Ja). The square 
phase is oriented such that the {10} planes are in contact with 
the liquid phase. Since each peak in the two particle correla- 
tion function represents a single plane, the surface energy of 
the interface is controlled by the width of the corresponding 
peak. The square phase {10}-facet has its surface energy de- 
rived from the width a-2 of the second peak of the correlation 
function, while the triangular phase facets are derived directly 
from the width «i of the first peak. The square phase sur- 
face energies can be tuned to isotropy by selecting the proper 
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(a) time = 3 



(b) time = 3.2 



(c) time = 3.5 



(d) time = 4 



FIG. 4. (Colour Online) 2D nucleation: from triangles to squares, 
(a) A liquid with density n — 0.09 is quenched into the triangular 
region of the phase diagram (a = 0.82), marked (A) in Fig.JUa), 
produces grains with triangular symmetry, (b-d) A second quench 
into the square portion of the phase diagram (er = 0.6), marked (B) in 
Fig.[5Ja), illustrating the nucleation and subsequent growth of square 
grains at the vertices of the triangular grains. 



ratio of the peak widths. Conversely, anisotropic surface en- 
ergy can be introduced by increasing or decreasing ol\ while 
holding oi2 constant. This effect also extends to the solid-solid 
surfaces between boundaries of different structure (for exam- 
ple, triangle-square boundaries) and misoriented solids of a 
structure derived from multiple peaks (for example, the square 
structure). The triangle structure is derived directly from a sin- 
gle peak and is therefore isotropic for a given misorientation. 
Small changes in the strength of the anisotropy lead to small 
changes in the angle of the solid liquid interfaces at the vertex 
point between the three phases. 

The second application of interest is the nucleation of a new 
structural phase by a quench into the second phase region of 
the phase diagram. For illustrative purposes, the simple 2D 
case of square-triangle structural transformations is used here. 
A liquid with a homogenized density field can be quenched 
into the triangular region of the phase diagram (see point A 
on the transformation path in Figure [3ja)) . Conserved noise 
in the dynamics is added, and seeds of the triangular phase 
nucleate and are allowed to coarsen for some time. One such 
small grain is illustrated in Fig. IS a). After some time, the 
material can be further quenched into the square portion of the 
phase diagram (see point B). Nucleation of the square phase 



and subsequent growth and coarsening of the new phase is 
observed. Note that the nucleation events happen at the vertex 
positions in the grain structures, as observed in experiments 
111]. Two such nucleation points are shown in Fig.Htb). The 
orientation of a nucleus and coherency strains can sometimes 
inhibit the growth of the new phase. This can cause the new 
phase to be dominated by another nucleation that is in a more 
favourable position to grow, an effect illustrated in panels (c) 
and (d) of Fig. S 

We have introduced new CDFT/PFC free energy function- 
al for efficient numerical study of solid state transformations. 
In contrast to earlier work [13J, we find that two-point correla- 
tions are sufficient to generate stable cubic lattices, and higher 
order terms in the expansion of AF ex can still be neglected. 
The correlation functions are systematically built up from fun- 
damental principles and desired crystallographic properties of 
phases of interest at a finite temperature. Our model captures 
a peritectic transition as well as the nucleation and growth of 
second-phase precipitates with different crystalline structures. 

This work has been supported by the Natural Science and 
Engineering Research Council of Canada (NSERC). 
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